## Load libraries


suppressPackageStartupMessages({
library(ArchR)
library(rhdf5)
library(tidyverse)
library(reticulate)
library(zellkonverter)
library(Matrix)
library(dichromat)
library(Seurat)
#library(caret)
h5disableFileLocking()})
rna_seurat <- readRDS("Seurat_objects/rna_Seurat_object")
hvg <- VariableFeatures(rna_seurat)

Prepare Gene expression matrix:

# get gene expression matrix
#proj <- loadArchRProject("12_Ricards_peaks_ChromVar")
proj <- loadArchRProject("/omics/groups/OE0533/internal/katharina/scDoRI/gastrulation_data/12_Ricards_peaks_ChromVar/")
metadata <- as.data.frame(getCellColData(proj))



gene_expr <- getMatrixFromProject(proj, useMatrix = "GeneExpressionMatrix")
gene_expr_mat <- assays(gene_expr)[[1]]
rownames(gene_expr_mat) <- rowData(gene_expr)$name
gene_expr_mat[1:5, 1:5]
## 5 x 5 sparse Matrix of class "dgCMatrix"
##        E8.5_rep1#TTACGTTTCTGGCATG-1 E8.5_rep1#TCCTCTAAGTCCTTCA-1
## Xkr4                      .                             .       
## Rp1                       .                             .       
## Sox17                     .                             .       
## Mrpl15                    0.8766163                     0.846525
## Lypla1                    0.8766163                     .       
##        E8.5_rep1#TATCCAGCACAGACTC-1 E8.5_rep1#GAGAACCAGACACTTA-1
## Xkr4                      .                                    .
## Rp1                       .                                    .
## Sox17                     0.2216017                            .
## Mrpl15                    0.6648052                            .
## Lypla1                    2.6592208                            .
##        E8.5_rep1#TCAAGTATCTTAGCGG-1
## Xkr4                        .      
## Rp1                         .      
## Sox17                       .      
## Mrpl15                      2.53893
## Lypla1                      .


#normalize gene expression
lognorm <- t(t(gene_expr_mat) / colSums(gene_expr_mat))
lognorm <- log(lognorm * 1e4 + 1)
#
# Gata1 <- lognorm[rownames(lognorm) %in% c("Gata1"), ]
# Sox9 <- lognorm[rownames(lognorm) %in% c("Sox9"), ]


# get the metadata


# metadata["Gata1"] <- Gata1
# metadata["Sox9"] <- Sox9
# metadata %>% head

Prepare motif deviation z-score matrix:

deviations <- getVarDeviations(proj, name = "MotifMatrix", plot = FALSE)
## DataFrame with 6 rows and 6 columns
##      seqnames     idx      name combinedVars combinedMeans      rank
##         <Rle> <array>   <array>    <numeric>     <numeric> <integer>
## f382        z     382 Gata6_382      61.0523      0.410885         1
## f386        z     386 Gata4_386      60.4598      0.407792         2
## f385        z     385 Gata5_385      59.1621      0.385362         3
## f387        z     387 Gata1_387      56.0319      0.365879         4
## f384        z     384 Gata3_384      52.9865      0.374011         5
## f383        z     383 Gata2_383      52.8551      0.366560         6
deviation_z_scoes <- deviations %>% as.data.frame() %>% filter(seqnames == "z") 

# get motif matrix
motifs <- getMatrixFromProject(proj, useMatrix = "MotifMatrix")
motif_mtx <- assays(motifs)[[2]]
# remove index number from TFs
tfs <- str_remove(rownames(motif_mtx), "_(?=[0-9])")
rownames(motif_mtx) <- tfs


# # get only one tf from matrix
# gata1_motif <- motif_mtx[rownames(motif_mtx) ==  tfs[grepl("^Gata1", tfs)], ]
# metadata["gata1_z_scores"] = gata1_motif
# 
# sox9_motif <- motif_mtx[rownames(motif_mtx) ==  tfs[grepl("^Sox9", tfs)], ]
# metadata["sox9_z_scores"] = sox9_motif

Prepare gene activity scores computed from p2g links:

# proj1 <- loadArchRProject("14_gene_scores_from_p2g_as_gene_expr_matrix/")
# gene_scores <- getMatrixFromProject(proj1, useMatrix = "GeneExpressionMatrix")
# gene_scores_mat <- assays(gene_scores)[[1]]
# rownames(gene_scores_mat) <- rowData(gene_scores)$name


gene_scores <- getMatrixFromProject(proj, useMatrix = "GeneScoreMatrix")
gene_scores_mat <- assays(gene_scores)[[1]]
rownames(gene_scores_mat) <- rowData(gene_scores)$name
colnames(gene_scores_mat) <- colnames(gene_scores)
# use new gene scores based on positive peak2gene links
#gene_scores <- readH5AD("jupyter_notebooks/p2g_gene_activity_scores/z_score_p2g_gene_activity_scores")
#gene_scores_mat <- assays(gene_scores)[[1]]






# add gene scores for gata1 to dataframe
# gata1_score_p2g <- gene_scores_mat[rownames(gene_scores_mat) %in% c("Gata1"), ]
# metadata["gata1_score_p2g"] <- gata1_score_p2g
# 
# # add gene scores for sox9 to dataframe
# sox9_score_p2g <- gene_scores_mat[rownames(gene_scores_mat) %in% c("Sox9"), ]
# metadata["sox9_score_p2g"] <- sox9_score_p2g

Prepare motif deviation scores computed using SEACells:

dev <- readRDS("/omics/groups/OE0533/internal/katharina/scDoRI/gastrulation_data/jupyter_notebooks/SEACell_files/SEA_aggregates_to_R/SEACell_ChromVarDev")


sea_archr_meta <- read_csv("/omics/groups/OE0533/internal/katharina/scDoRI/gastrulation_data/jupyter_notebooks/SEACell_files/archr_sea_metadata.csv")

celltype_df <- sea_archr_meta %>%
  dplyr::count(SEACell, celltypes) %>%
  dplyr::group_by(SEACell) %>%
  slice_max(order_by=n, n=1, with_ties = FALSE) %>% 
  select(SEACell, celltypes)


df <- colData(dev) %>% as.data.frame() %>% rownames_to_column("SEACell")
df <- left_join(celltype_df, df, by = "SEACell")


colData(dev) <- DataFrame(df)

sea_mtx <- assays(dev)[[1]]


sea_meta <- colData(dev) %>% as.data.frame()
#tfs <- rownames(dev)
#metadata <- colData(dev) %>% as.data.frame()
proj1 <- loadArchRProject("/omics/groups/OE0533/internal/katharina/scDoRI/gastrulation_data_new/Kathi/06_deep_chromvar/")
# get motif matrix
deep_motifs <- getMatrixFromProject(proj1, useMatrix = "DeepLearningMotifs1")
deep_motif_mtx <- assays(deep_motifs)[[2]]

rownames(deep_motif_mtx) <- tolower(rownames(deep_motif_mtx))
substr(rownames(deep_motif_mtx), 1, 1) <- toupper(substr(rownames(deep_motif_mtx), 1, 1))

Color Palette:

colPalette_celltypes = c('#532C8A',
 '#c19f70',
 '#f9decf',
 '#c9a997',
 '#B51D8D',
 '#3F84AA',
 '#9e6762',
 '#354E23',
 '#F397C0',
 '#ff891c',
 '#635547',
 '#C72228',
 '#f79083',
 '#EF4E22',
 '#989898',
 '#7F6874',
 '#8870ad',
 '#647a4f',
 '#EF5A9D',
 '#FBBE92',
 '#139992',
 '#cc7818',
 '#DFCDE4',
 '#8EC792',
 '#C594BF',
 '#C3C388',
 '#0F4A9C',
 '#FACB12',
 '#8DB5CE',
 '#1A1A1A',
 '#C9EBFB',
 '#DABE99',
 '#65A83E',
 '#005579',
 '#CDE088',
 '#f7f79e',
 '#F6BFCB')

celltypes <- (as.data.frame(getCellColData(proj)) %>% group_by(celltypes) %>% 
 summarise(n = n()))$celltypes

col <- setNames(colPalette_celltypes, celltypes)

Select a few marker genes and transcription factors:

marker_genes <- c("Lamb1",  "Sparc", "Elf5", "Ascl2", "Tfap2c", "Ttr",
                  "Apoa2", "Apoe", "Cystm1", "Emb", "Spink1",  "Krt19",
                  "Dkk1", "Grhl3", "Trp63", "Grhl2",  "Pax6", "Pax2",
                  "En1", "Foxd3", "Tfap2a", "Pax3", "Sox9",
                  "Six3", "Hesx1", "Irx3", "Hoxb9", "Cdx4",
                  "Hes3", "Hba-a2", "Hba-a1",  "Hbb-bh1", "Gata1",
                   "Gata6",
                  #"Gata6", "Gata5",
                  "Cited4",
                   "Cdh5", "Pecam1", "Anxa5", "Etv2", "Igf2",
                  "Krt8", "Krt18", "Pmp22", "Ahnak", "Bmp4", "Tbx4", "Hoxa11",
                  "Hoxa10", "Tnnt2", "Myl4",  "Myl7", "Acta2",
                  "Smarcd3", "Tcf21", "Tbx6", "Dll1", "Aldh1a2", "Tcf15",
                  "Meox1", "Tbx1", "Gbx2", "Cdx1", "Hoxb1", "Hes7", "Osr1",
                  "Mesp2", "Lefty2", "Mesp1", "Cer1",  "Chrd", 
                  "Foxa2", "Pax7", "Fgf8", "Lhx1", "Mixl1", "Otx2", "Hhex",
                   "Ifitm3", "Nkx1-2", "Eomes", "Nanog", "Utf1",
                  "Epcam", "Pou5f1" )
#"Sox2"
#"Gata2"
#"Gata4"
#  "Gata5",
#"Gata6",
# "Gsc",
  p <- metadata %>%
    mutate(!!n := motif_n) %>%
    group_by(celltypes) %>%
    #summarise_at(vars(n), funs(mean(., na.rm=TRUE)))
    summarise(mean = mean(!!(sym(n)))) %>%
    ggplot() +
    geom_bar(aes(x = celltypes, y = mean, fill = celltypes), stat = "identity") +
    scale_fill_manual(values = col) +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), legend.position = "None") +#%>% print()
    labs(title = paste0(n), y = "SEACell deviation score")
  print(p)

GATA Factors

for (n in c("Gata1", "Gata2", "Gata3", "Gata4", "Gata6")) { # "Gata5",
  print(n)
  # select one row for a particular gene
  score_n <- gene_scores_mat[rownames(gene_scores_mat) %in%  n, ]
  # add score for this gene to metadata
  metadata[paste0("score_",n)] <- score_n
  # select gene expression for a particular gene
  #expr_n <- lognorm[rownames(lognorm) %in% c(n), ]
  #metadata[paste0("expr_", n)] <- expr_n
  
  
  seacells_n <- deep_motif_mtx[n, ]
  metadata[paste0("seacell_", n)] <- seacells_n

  # select motif z score
  motif_n <- motif_mtx[rownames(motif_mtx) ==  tfs[grepl(paste0("^", n), tfs)], ]
  metadata[paste0("motif_", n)] <- motif_n
  
  # create barplots for gene scores, gene expression and motif z-score
  plots <- map(c("score_", "motif_"), function(p){
    df <- metadata %>%
    group_by(celltypes) %>%
    summarise_at(vars(paste0(p, n)), funs(mean(., na.rm=TRUE)))
    df %>% ggplot() +
    geom_bar(aes(x = celltypes, y = df %>% pull(paste0(p, n)), 
                 fill = celltypes), stat = "identity") +
    scale_fill_manual(values = col) +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
          legend.position = "none") +
    labs(title = paste0(n), x = "celltype", y = paste0(p))
  })
  sea_plot <- map(seq.int(1), function(sea){
    metadata %>% 
    group_by(celltypes) %>%
    summarise(mean = mean(!!sym(paste0("seacell_",n)))) %>%
    ggplot() +
    geom_bar(aes(x = celltypes, y = mean, fill = celltypes), stat = "identity") +
    scale_fill_manual(values = col) +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), legend.position = "None") +#%>% print()
    labs(title = paste0(n), y = "SEACell deviation score")
  })
  # create scatter plots for gene expression and motif z-score
  scatter_plots <- map(seq.int(1), function(s){
    df <- metadata %>% group_by(celltypes) %>% 
    summarise_at(vars(paste0("motif_", n), 
                      paste0("score_", n)), 
                 funs(mean(., na.rm = TRUE)))
    df %>%
    ggplot() +
    geom_smooth(aes(x = df %>% pull(paste0("score_", n)), 
                    y = df %>% pull(paste0("motif_", n))),
                formula = y ~ x, method = "lm", size = .1) +
    geom_point(aes(x = df %>% pull(paste0("score_", n)), 
                   y = df %>% pull(paste0("motif_", n)), 
                   col = celltypes, size = 1)) +
    #labs(x = "Gata1 gene expression", y = "Gata1 motif accessibility (z-score)") +
    theme(legend.position = "none") +
    scale_color_manual(values = col) +
    labs(title = paste0(n), x = "gene activity score", y = "TF-motif z-score")
    })
  do.call(gridExtra::grid.arrange, c(plots, sea_plot, scatter_plots, ncol = 2))

}
## [1] "Gata1"
## Warning: `funs()` was deprecated in dplyr 0.8.0.
## Please use a list of either functions or lambdas: 
## 
##   # Simple named list: 
##   list(mean = mean, median = median)
## 
##   # Auto named with `tibble::lst()`: 
##   tibble::lst(mean, median)
## 
##   # Using lambdas
##   list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was generated.

## [1] "Gata2"

## [1] "Gata3"

## [1] "Gata4"

## [1] "Gata6"

Marker Genes

SEACells


for (n in marker_genes) {
  print(n)
  # select one row for a particular gene
  score_n <- gene_scores_mat[rownames(gene_scores_mat) %in% c(n), ]
  # add score for this gene to metadata
  metadata[paste0("score_",n)] <- score_n
  # select gene expression for a particular gene
  expr_n <- lognorm[rownames(lognorm) %in% c(n), ]
  metadata[paste0("expr_", n)] <- expr_n
  
  
  seacells_n <- sea_mtx[rownames(sea_mtx) == tfs[grepl(paste0("^", n), tfs)], ]
  sea_meta[paste0("seacell_", n)] <- seacells_n

  # if the marker gene is a TF
  if (length(tfs[grepl(paste0("^", n), tfs)]) > 0) {
    
      # select motif z score
      motif_n <- motif_mtx[rownames(motif_mtx) ==  tfs[grepl(paste0("^", n), tfs)], ]
      metadata[paste0("motif_", n)] <- motif_n
      
      # create barplots for gene scores, gene expression and motif z-score
      plots <- map(c("score_", "motif_"), function(p){
        df <- metadata %>%
        group_by(celltypes) %>%
        summarise_at(vars(paste0(p, n)), funs(mean(., na.rm=TRUE)))
        df %>% ggplot() +
        geom_bar(aes(x = celltypes, y = df %>% pull(paste0(p, n)), 
                     fill = celltypes), stat = "identity") +
        scale_fill_manual(values = col) +
        theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
              legend.position = "none") +
        labs(title = paste0(n), x = "celltype", y = paste0(p))
      })
      
      # create scatter plots for gene expression and motif z-score
      scatter_plots <- map(seq.int(1), function(s){
        df <- metadata %>% group_by(celltypes) %>% 
        summarise_at(vars(paste0("motif_", n), 
                          paste0("score_", n)), 
                     funs(mean(., na.rm = TRUE)))
        df %>%
        ggplot() +
        geom_smooth(aes(x = df %>% pull(paste0("score_", n)), 
                        y = df %>% pull(paste0("motif_", n))),
                    formula = y ~ x, method = "lm", size = .1) +
        geom_point(aes(x = df %>% pull(paste0("score_", n)), 
                       y = df %>% pull(paste0("motif_", n)), 
                       col = celltypes, size = 1)) +
        #labs(x = "Gata1 gene expression", y = "Gata1 motif accessibility (z-score)") +
        theme(legend.position = "none") +
        scale_color_manual(values = col) +
        labs(title = paste0(n), x = "gene activity score", y = "TF-motif z-score")
        })
        sea_plot <- map(seq.int(1), function(sea){
      sea_meta %>% 
      group_by(celltypes) %>%
      summarise(mean = mean(!!sym(paste0("seacell_",n)))) %>%
      ggplot() +
      geom_bar(aes(x = celltypes, y = mean, fill = celltypes), stat = "identity") +
      scale_fill_manual(values = col) +
      theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), legend.position = "None") +#%>% print()
      labs(title = paste0(n), y = "SEACell deviation score")
    })
      
    # if the marker gene is no TF
    } else {
      
      # create barplots for gene scores, gene expression and motif z-score
      plots <- map(c("score_"), function(p){
        df <- metadata %>%
        group_by(celltypes) %>%
        summarise_at(vars(paste0(p, n)), funs(mean(., na.rm=TRUE)))
        df %>% ggplot() +
        geom_bar(aes(x = celltypes, y = df %>% pull(paste0(p, n)), 
                     fill = celltypes), stat = "identity") +
        scale_fill_manual(values = col) +
        theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
              legend.position = "none") +
        labs(title = paste0(n), x = "celltype", y = paste0(p))

      })
      
      # create scatter plots for gene expression and gene_score
      scatter_plots <- map(seq.int(1), function(s){
        df <- metadata %>% group_by(celltypes) %>%
        summarise_at(vars(paste0("score_", n),
                          paste0("expr_", n)),
                     funs(mean(., na.rm = TRUE))) 
        df %>%
        ggplot() +
        geom_smooth(aes(x = df %>% pull(paste0("expr_", n)), 
                        y = df %>% pull(paste0("score_", n))),
                    formula = y ~ x, method = "lm", size = .1) +
        geom_point(aes(x = df %>% pull(paste0("expr_", n)),
                       y = df %>% pull(paste0("score_", n)),
                       col = celltypes, size = 1)) +
        #labs(x = "Gata1 gene expression", y = "Gata1 motif accessibility (z-score)") +
        theme(legend.position = "none") +
        scale_color_manual(values = col)+
        labs(title = paste0(n), x = "gene expression", y = "gene activity score")

      })
    }

  do.call(gridExtra::grid.arrange, c(plots, scatter_plots, ncol = 2, nrow = 2))
  # 

}
## [1] "Lamb1"

## [1] "Sparc"

## [1] "Elf5"

## [1] "Ascl2"

## [1] "Tfap2c"

## [1] "Ttr"

## [1] "Apoa2"

## [1] "Apoe"

## [1] "Cystm1"

## [1] "Emb"

## [1] "Spink1"

## [1] "Krt19"

## [1] "Dkk1"

## [1] "Grhl3"

## [1] "Trp63"

## [1] "Grhl2"

## [1] "Pax6"

## [1] "Pax2"

## [1] "En1"

## [1] "Foxd3"

## [1] "Tfap2a"

## [1] "Pax3"

## [1] "Sox9"

## [1] "Six3"

## [1] "Hesx1"

## [1] "Irx3"

## [1] "Hoxb9"

## [1] "Cdx4"

## [1] "Hes3"

## [1] "Hba-a2"

## [1] "Hba-a1"

## [1] "Hbb-bh1"

## [1] "Gata1"

## [1] "Gata6"

## [1] "Cited4"

## [1] "Cdh5"

## [1] "Pecam1"

## [1] "Anxa5"

## [1] "Etv2"

## [1] "Igf2"

## [1] "Krt8"

## [1] "Krt18"

## [1] "Pmp22"

## [1] "Ahnak"

## [1] "Bmp4"

## [1] "Tbx4"

## [1] "Hoxa11"

## [1] "Hoxa10"

## [1] "Tnnt2"

## [1] "Myl4"

## [1] "Myl7"

## [1] "Acta2"

## [1] "Smarcd3"

## [1] "Tcf21"

## [1] "Tbx6"

## [1] "Dll1"

## [1] "Aldh1a2"

## [1] "Tcf15"

## [1] "Meox1"
## [1] "Tbx1"
## Warning in rownames(sea_mtx) == tfs[grepl(paste0("^", n), tfs)]: longer object
## length is not a multiple of shorter object length
## Warning in rownames(motif_mtx) == tfs[grepl(paste0("^", n), tfs)]: longer object
## length is not a multiple of shorter object length

## [1] "Gbx2"

## [1] "Cdx1"

## [1] "Hoxb1"

## [1] "Hes7"

## [1] "Osr1"

## [1] "Mesp2"

## [1] "Lefty2"

## [1] "Mesp1"

## [1] "Cer1"

## [1] "Chrd"

## [1] "Foxa2"

## [1] "Pax7"

## [1] "Fgf8"

## [1] "Lhx1"

## [1] "Mixl1"

## [1] "Otx2"

## [1] "Hhex"

## [1] "Ifitm3"

## [1] "Nkx1-2"

## [1] "Eomes"

## [1] "Nanog"

## [1] "Utf1"

## [1] "Epcam"

## [1] "Pou5f1"

Plot gene expression, gene scores & motif z-scores:

```#{r, fig.width=10, fig.height=10}

for (n in marker_genes) { print(n) # select one row for a particular gene score_n <- gene_scores_mat[rownames(gene_scores_mat) %in% c(n), ] # add score for this gene to metadata metadata[paste0(“score_”,n)] <- score_n # select gene expression for a particular gene expr_n <- lognorm[rownames(lognorm) %in% c(n), ] metadata[paste0(“expr_”, n)] <- expr_n

seacells_n <- sea_mtx[rownames(sea_mtx) == tfs[grepl(paste0(“^”, n), tfs)], ] sea_meta[paste0(“seacell_”, n)] <- seacells_n

# if the marker gene is a TF if (length(tfs[grepl(paste0(“^”, n), tfs)]) > 0) {

  # select motif z score
  motif_n <- motif_mtx[rownames(motif_mtx) ==  tfs[grepl(paste0("^", n), tfs)], ]
  metadata[paste0("motif_", n)] <- motif_n
  
  # create barplots for gene scores, gene expression and motif z-score
  plots <- map(c("score_", "expr_", "motif_"), function(p){
    df <- metadata %>%
    group_by(celltypes) %>%
    summarise_at(vars(paste0(p, n)), funs(mean(., na.rm=TRUE)))
    df %>% ggplot() +
    geom_bar(aes(x = celltypes, y = df %>% pull(paste0(p, n)), 
                 fill = celltypes), stat = "identity") +
    scale_fill_manual(values = col) +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
          legend.position = "none") +
    labs(title = paste0(n), x = "celltype", y = paste0(p))
  })
  
  # create scatter plots for gene expression and motif z-score
  scatter_plots <- map(seq.int(1), function(s){
    df <- metadata %>% group_by(celltypes) %>% 
    summarise_at(vars(paste0("motif_", n), 
                      paste0("score_", n)), 
                 funs(mean(., na.rm = TRUE)))
    df %>%
    ggplot() +
    geom_smooth(aes(x = df %>% pull(paste0("score_", n)), 
                    y = df %>% pull(paste0("motif_", n))),
                formula = y ~ x, method = "lm", size = .1) +
    geom_point(aes(x = df %>% pull(paste0("score_", n)), 
                   y = df %>% pull(paste0("motif_", n)), 
                   col = celltypes, size = 1)) +
    #labs(x = "Gata1 gene expression", y = "Gata1 motif accessibility (z-score)") +
    theme(legend.position = "none") +
    scale_color_manual(values = col) +
    labs(title = paste0(n), x = "gene activity score", y = "TF-motif z-score")
    })
  
# if the marker gene is no TF
} else {
  
  # create barplots for gene scores, gene expression and motif z-score
  plots <- map(c("score_", "expr_"), function(p){
    df <- metadata %>%
    group_by(celltypes) %>%
    summarise_at(vars(paste0(p, n)), funs(mean(., na.rm=TRUE)))
    df %>% ggplot() +
    geom_bar(aes(x = celltypes, y = df %>% pull(paste0(p, n)), 
                 fill = celltypes), stat = "identity") +
    scale_fill_manual(values = col) +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
          legend.position = "none") +
    labs(title = paste0(n), x = "celltype", y = paste0(p))

  })
  
  # create scatter plots for gene expression and gene_score
  scatter_plots <- map(seq.int(1), function(s){
    df <- metadata %>% group_by(celltypes) %>%
    summarise_at(vars(paste0("score_", n),
                      paste0("expr_", n)),
                 funs(mean(., na.rm = TRUE))) 
    df %>%
    ggplot() +
    geom_smooth(aes(x = df %>% pull(paste0("expr_", n)), 
                    y = df %>% pull(paste0("score_", n))),
                formula = y ~ x, method = "lm", size = .1) +
    geom_point(aes(x = df %>% pull(paste0("expr_", n)),
                   y = df %>% pull(paste0("score_", n)),
                   col = celltypes, size = 1)) +
    #labs(x = "Gata1 gene expression", y = "Gata1 motif accessibility (z-score)") +
    theme(legend.position = "none") +
    scale_color_manual(values = col)+
    labs(title = paste0(n), x = "gene expression", y = "gene activity score")

  })
}

do.call(gridExtra::grid.arrange, c(plots, scatter_plots, ncol = 2, nrow = 2)) #

}




```#{r, fig.width=10, fig.height=10}
p1 <- metadata %>%
  group_by(celltypes) %>%
  summarise_at(vars(Gata1), funs(median(., na.rm=TRUE))) %>% ggplot() +
  geom_bar(aes(x = celltypes, y = Gata1, fill = celltypes), stat = "identity") +
  scale_fill_manual(values = col) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), 
        legend.position = "none")

p2 <- metadata %>%
  group_by(celltypes) %>%
  summarise_at(vars(gata1_score_p2g), funs(mean(., na.rm=TRUE))) %>% ggplot() +
  geom_bar(aes(x = celltypes, y = gata1_score_p2g, fill = celltypes), stat = "identity") +
  scale_fill_manual(values = col) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), 
        legend.position = "none")


p3 <- metadata %>% group_by(celltypes) %>% 
  summarise_at(vars(gata1_z_scores, Gata1), funs(mean(., na.rm = TRUE))) %>%
  ggplot() +
  geom_smooth(aes(x = Gata1, y = gata1_z_scores), formula = y ~ x, method = "lm", size = .1) +
  geom_point(aes(x = Gata1, y = gata1_z_scores, col = celltypes, size = 1)) +
  labs(x = "Gata1 gene expression", y = "Gata1 motif accessibility (z-score)") +
  theme(legend.position = "none") +
  scale_color_manual(values = col) 

p4 <- metadata %>%
  group_by(celltypes) %>%
  summarise_at(vars(gata1_z_scores), funs(mean(., na.rm=TRUE))) %>% ggplot() +
  geom_bar(aes(x = celltypes, y = gata1_z_scores, fill = celltypes), stat = "identity") +
  scale_fill_manual(values = col) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), 
        legend.position = "none")

cowplot::plot_grid(p1, p2, p3, p4)

p1 <- metadata %>%
  group_by(celltypes) %>%
  summarise_at(vars(Sox9), funs(mean(., na.rm=TRUE))) %>% ggplot() +
  geom_bar(aes(x = celltypes, y = Sox9, fill = celltypes), stat = "identity") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1),
        legend.position = "none") +
  scale_fill_manual(values = col)



p2 <- metadata %>%
  group_by(celltypes) %>%
  summarise_at(vars(sox9_z_scores), funs(mean(., na.rm=TRUE))) %>% ggplot() +
  geom_bar(aes(x = celltypes, y = sox9_z_scores, fill = celltypes), stat = "identity") +
  scale_fill_manual(values = col) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), 
        legend.position = "none")

p3 <- metadata %>%
  group_by(celltypes) %>%
  summarise_at(vars(sox9_score_p2g), funs(mean(., na.rm=TRUE))) %>% ggplot() +
  geom_bar(aes(x = celltypes, y = sox9_score_p2g, fill = celltypes), stat = "identity") +
  scale_fill_manual(values = col) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), 
        legend.position = "none")


p4 <- metadata %>% group_by(celltypes) %>% 
  summarise_at(vars(sox9_z_scores, Sox9), funs(mean(., na.rm = TRUE))) %>%
  ggplot() +
  geom_smooth(aes(x = Sox9, y = sox9_z_scores), formula = y ~ x, method = "lm", size = .1) +
  geom_point(aes(x = Sox9, y = sox9_z_scores, col = celltypes, size = 1)) +
  labs(x = "Gata1 gene expression", y = "Gata1 motif accessibility (z-score)") +
  theme(legend.position = "none") +
  scale_color_manual(values = col) 

cowplot::plot_grid(p1, p2, p3, p4)
LS0tCnRpdGxlOiAiMTRfQmFycGxvdHMiCm91dHB1dDogCiAgaHRtbF9kb2N1bWVudDoKICAgIHRvYzogdHJ1ZQogICAgdG9jX2RlcHRoOiA1CiAgICBjb2RlX2ZvbGRpbmc6IGhpZGUKICAgIHRvY19mbG9hdDogdHJ1ZQogICAgY29kZV9kb3dubG9hZDogdHJ1ZQogICAgdGhlbWU6IGNvc21vCiAgICBoaWdobGlnaHQ6IHRleHRtYXRlCi0tLQoKPHN0eWxlPgpib2R5IHsKdGV4dC1hbGlnbjoganVzdGlmeX0KPC9zdHlsZT4KCmBgYHtyIHNldHVwLCBpbmNsdWRlPUZBTFNFfQprbml0cjo6b3B0c19jaHVuayRzZXQoY2FjaGUgPSBGQUxTRSwgYXV0b2RlcCA9IFRSVUUsIAogICAgICAgICAgICAgICAgICAgICAgY29sbGFwc2UgPSBUUlVFLCBtZXNzYWdlID0gRkFMU0UpCmtuaXRyOjpvcHRzX2tuaXQkc2V0KHJvb3QuZGlyID0gIi9vbWljcy9ncm91cHMvT0UwNTMzL2ludGVybmFsL2thdGhhcmluYS9zY0RvUkkvZ2FzdHJ1bGF0aW9uX2RhdGEvIikKc2V0d2QoIi9vbWljcy9ncm91cHMvT0UwNTMzL2ludGVybmFsL2thdGhhcmluYS9zY0RvUkkvZ2FzdHJ1bGF0aW9uX2RhdGEvIikKCnNldC5zZWVkKDEpCmBgYAoKICAKYGBge3J9CiMjIExvYWQgbGlicmFyaWVzCgoKc3VwcHJlc3NQYWNrYWdlU3RhcnR1cE1lc3NhZ2VzKHsKbGlicmFyeShBcmNoUikKbGlicmFyeShyaGRmNSkKbGlicmFyeSh0aWR5dmVyc2UpCmxpYnJhcnkocmV0aWN1bGF0ZSkKbGlicmFyeSh6ZWxsa29udmVydGVyKQpsaWJyYXJ5KE1hdHJpeCkKbGlicmFyeShkaWNocm9tYXQpCmxpYnJhcnkoU2V1cmF0KQojbGlicmFyeShjYXJldCkKaDVkaXNhYmxlRmlsZUxvY2tpbmcoKX0pCmBgYAoKYGBgI3tyfQpybmFfc2V1cmF0IDwtIHJlYWRSRFMoIlNldXJhdF9vYmplY3RzL3JuYV9TZXVyYXRfb2JqZWN0IikKaHZnIDwtIFZhcmlhYmxlRmVhdHVyZXMocm5hX3NldXJhdCkKYGBgCgpQcmVwYXJlIEdlbmUgZXhwcmVzc2lvbiBtYXRyaXg6CgpgYGB7cn0KIyBnZXQgZ2VuZSBleHByZXNzaW9uIG1hdHJpeAojcHJvaiA8LSBsb2FkQXJjaFJQcm9qZWN0KCIxMl9SaWNhcmRzX3BlYWtzX0Nocm9tVmFyIikKcHJvaiA8LSBsb2FkQXJjaFJQcm9qZWN0KCIvb21pY3MvZ3JvdXBzL09FMDUzMy9pbnRlcm5hbC9rYXRoYXJpbmEvc2NEb1JJL2dhc3RydWxhdGlvbl9kYXRhLzEyX1JpY2FyZHNfcGVha3NfQ2hyb21WYXIvIikKbWV0YWRhdGEgPC0gYXMuZGF0YS5mcmFtZShnZXRDZWxsQ29sRGF0YShwcm9qKSkKCgoKZ2VuZV9leHByIDwtIGdldE1hdHJpeEZyb21Qcm9qZWN0KHByb2osIHVzZU1hdHJpeCA9ICJHZW5lRXhwcmVzc2lvbk1hdHJpeCIpCmdlbmVfZXhwcl9tYXQgPC0gYXNzYXlzKGdlbmVfZXhwcilbWzFdXQpyb3duYW1lcyhnZW5lX2V4cHJfbWF0KSA8LSByb3dEYXRhKGdlbmVfZXhwcikkbmFtZQpnZW5lX2V4cHJfbWF0WzE6NSwgMTo1XQoKCiNub3JtYWxpemUgZ2VuZSBleHByZXNzaW9uCmxvZ25vcm0gPC0gdCh0KGdlbmVfZXhwcl9tYXQpIC8gY29sU3VtcyhnZW5lX2V4cHJfbWF0KSkKbG9nbm9ybSA8LSBsb2cobG9nbm9ybSAqIDFlNCArIDEpCiMKIyBHYXRhMSA8LSBsb2dub3JtW3Jvd25hbWVzKGxvZ25vcm0pICVpbiUgYygiR2F0YTEiKSwgXQojIFNveDkgPC0gbG9nbm9ybVtyb3duYW1lcyhsb2dub3JtKSAlaW4lIGMoIlNveDkiKSwgXQoKCiMgZ2V0IHRoZSBtZXRhZGF0YQoKCiMgbWV0YWRhdGFbIkdhdGExIl0gPC0gR2F0YTEKIyBtZXRhZGF0YVsiU294OSJdIDwtIFNveDkKIyBtZXRhZGF0YSAlPiUgaGVhZAoKCgpgYGAKClByZXBhcmUgbW90aWYgZGV2aWF0aW9uIHotc2NvcmUgbWF0cml4OgoKYGBge3J9CmRldmlhdGlvbnMgPC0gZ2V0VmFyRGV2aWF0aW9ucyhwcm9qLCBuYW1lID0gIk1vdGlmTWF0cml4IiwgcGxvdCA9IEZBTFNFKQpkZXZpYXRpb25fel9zY29lcyA8LSBkZXZpYXRpb25zICU+JSBhcy5kYXRhLmZyYW1lKCkgJT4lIGZpbHRlcihzZXFuYW1lcyA9PSAieiIpIAoKIyBnZXQgbW90aWYgbWF0cml4Cm1vdGlmcyA8LSBnZXRNYXRyaXhGcm9tUHJvamVjdChwcm9qLCB1c2VNYXRyaXggPSAiTW90aWZNYXRyaXgiKQptb3RpZl9tdHggPC0gYXNzYXlzKG1vdGlmcylbWzJdXQojIHJlbW92ZSBpbmRleCBudW1iZXIgZnJvbSBURnMKdGZzIDwtIHN0cl9yZW1vdmUocm93bmFtZXMobW90aWZfbXR4KSwgIl8oPz1bMC05XSkiKQpyb3duYW1lcyhtb3RpZl9tdHgpIDwtIHRmcwoKCiMgIyBnZXQgb25seSBvbmUgdGYgZnJvbSBtYXRyaXgKIyBnYXRhMV9tb3RpZiA8LSBtb3RpZl9tdHhbcm93bmFtZXMobW90aWZfbXR4KSA9PSAgdGZzW2dyZXBsKCJeR2F0YTEiLCB0ZnMpXSwgXQojIG1ldGFkYXRhWyJnYXRhMV96X3Njb3JlcyJdID0gZ2F0YTFfbW90aWYKIyAKIyBzb3g5X21vdGlmIDwtIG1vdGlmX210eFtyb3duYW1lcyhtb3RpZl9tdHgpID09ICB0ZnNbZ3JlcGwoIl5Tb3g5IiwgdGZzKV0sIF0KIyBtZXRhZGF0YVsic294OV96X3Njb3JlcyJdID0gc294OV9tb3RpZgpgYGAKCgpQcmVwYXJlIGdlbmUgYWN0aXZpdHkgc2NvcmVzIGNvbXB1dGVkIGZyb20gcDJnIGxpbmtzOgoKYGBge3J9CiMgcHJvajEgPC0gbG9hZEFyY2hSUHJvamVjdCgiMTRfZ2VuZV9zY29yZXNfZnJvbV9wMmdfYXNfZ2VuZV9leHByX21hdHJpeC8iKQojIGdlbmVfc2NvcmVzIDwtIGdldE1hdHJpeEZyb21Qcm9qZWN0KHByb2oxLCB1c2VNYXRyaXggPSAiR2VuZUV4cHJlc3Npb25NYXRyaXgiKQojIGdlbmVfc2NvcmVzX21hdCA8LSBhc3NheXMoZ2VuZV9zY29yZXMpW1sxXV0KIyByb3duYW1lcyhnZW5lX3Njb3Jlc19tYXQpIDwtIHJvd0RhdGEoZ2VuZV9zY29yZXMpJG5hbWUKCgpnZW5lX3Njb3JlcyA8LSBnZXRNYXRyaXhGcm9tUHJvamVjdChwcm9qLCB1c2VNYXRyaXggPSAiR2VuZVNjb3JlTWF0cml4IikKZ2VuZV9zY29yZXNfbWF0IDwtIGFzc2F5cyhnZW5lX3Njb3JlcylbWzFdXQpyb3duYW1lcyhnZW5lX3Njb3Jlc19tYXQpIDwtIHJvd0RhdGEoZ2VuZV9zY29yZXMpJG5hbWUKY29sbmFtZXMoZ2VuZV9zY29yZXNfbWF0KSA8LSBjb2xuYW1lcyhnZW5lX3Njb3JlcykKIyB1c2UgbmV3IGdlbmUgc2NvcmVzIGJhc2VkIG9uIHBvc2l0aXZlIHBlYWsyZ2VuZSBsaW5rcwojZ2VuZV9zY29yZXMgPC0gcmVhZEg1QUQoImp1cHl0ZXJfbm90ZWJvb2tzL3AyZ19nZW5lX2FjdGl2aXR5X3Njb3Jlcy96X3Njb3JlX3AyZ19nZW5lX2FjdGl2aXR5X3Njb3JlcyIpCiNnZW5lX3Njb3Jlc19tYXQgPC0gYXNzYXlzKGdlbmVfc2NvcmVzKVtbMV1dCgoKCgoKCiMgYWRkIGdlbmUgc2NvcmVzIGZvciBnYXRhMSB0byBkYXRhZnJhbWUKIyBnYXRhMV9zY29yZV9wMmcgPC0gZ2VuZV9zY29yZXNfbWF0W3Jvd25hbWVzKGdlbmVfc2NvcmVzX21hdCkgJWluJSBjKCJHYXRhMSIpLCBdCiMgbWV0YWRhdGFbImdhdGExX3Njb3JlX3AyZyJdIDwtIGdhdGExX3Njb3JlX3AyZwojIAojICMgYWRkIGdlbmUgc2NvcmVzIGZvciBzb3g5IHRvIGRhdGFmcmFtZQojIHNveDlfc2NvcmVfcDJnIDwtIGdlbmVfc2NvcmVzX21hdFtyb3duYW1lcyhnZW5lX3Njb3Jlc19tYXQpICVpbiUgYygiU294OSIpLCBdCiMgbWV0YWRhdGFbInNveDlfc2NvcmVfcDJnIl0gPC0gc294OV9zY29yZV9wMmcKYGBgCgpQcmVwYXJlIG1vdGlmIGRldmlhdGlvbiBzY29yZXMgY29tcHV0ZWQgdXNpbmcgU0VBQ2VsbHM6CgoKYGBge3J9CmRldiA8LSByZWFkUkRTKCIvb21pY3MvZ3JvdXBzL09FMDUzMy9pbnRlcm5hbC9rYXRoYXJpbmEvc2NEb1JJL2dhc3RydWxhdGlvbl9kYXRhL2p1cHl0ZXJfbm90ZWJvb2tzL1NFQUNlbGxfZmlsZXMvU0VBX2FnZ3JlZ2F0ZXNfdG9fUi9TRUFDZWxsX0Nocm9tVmFyRGV2IikKCgpzZWFfYXJjaHJfbWV0YSA8LSByZWFkX2NzdigiL29taWNzL2dyb3Vwcy9PRTA1MzMvaW50ZXJuYWwva2F0aGFyaW5hL3NjRG9SSS9nYXN0cnVsYXRpb25fZGF0YS9qdXB5dGVyX25vdGVib29rcy9TRUFDZWxsX2ZpbGVzL2FyY2hyX3NlYV9tZXRhZGF0YS5jc3YiKQoKY2VsbHR5cGVfZGYgPC0gc2VhX2FyY2hyX21ldGEgJT4lCiAgZHBseXI6OmNvdW50KFNFQUNlbGwsIGNlbGx0eXBlcykgJT4lCiAgZHBseXI6Omdyb3VwX2J5KFNFQUNlbGwpICU+JQogIHNsaWNlX21heChvcmRlcl9ieT1uLCBuPTEsIHdpdGhfdGllcyA9IEZBTFNFKSAlPiUgCiAgc2VsZWN0KFNFQUNlbGwsIGNlbGx0eXBlcykKCgpkZiA8LSBjb2xEYXRhKGRldikgJT4lIGFzLmRhdGEuZnJhbWUoKSAlPiUgcm93bmFtZXNfdG9fY29sdW1uKCJTRUFDZWxsIikKZGYgPC0gbGVmdF9qb2luKGNlbGx0eXBlX2RmLCBkZiwgYnkgPSAiU0VBQ2VsbCIpCgoKY29sRGF0YShkZXYpIDwtIERhdGFGcmFtZShkZikKCnNlYV9tdHggPC0gYXNzYXlzKGRldilbWzFdXQoKCnNlYV9tZXRhIDwtIGNvbERhdGEoZGV2KSAlPiUgYXMuZGF0YS5mcmFtZSgpCiN0ZnMgPC0gcm93bmFtZXMoZGV2KQojbWV0YWRhdGEgPC0gY29sRGF0YShkZXYpICU+JSBhcy5kYXRhLmZyYW1lKCkKCgpgYGAKCmBgYHtyfQpwcm9qMSA8LSBsb2FkQXJjaFJQcm9qZWN0KCIvb21pY3MvZ3JvdXBzL09FMDUzMy9pbnRlcm5hbC9rYXRoYXJpbmEvc2NEb1JJL2dhc3RydWxhdGlvbl9kYXRhX25ldy9LYXRoaS8wNl9kZWVwX2Nocm9tdmFyLyIpCiMgZ2V0IG1vdGlmIG1hdHJpeApkZWVwX21vdGlmcyA8LSBnZXRNYXRyaXhGcm9tUHJvamVjdChwcm9qMSwgdXNlTWF0cml4ID0gIkRlZXBMZWFybmluZ01vdGlmczEiKQpkZWVwX21vdGlmX210eCA8LSBhc3NheXMoZGVlcF9tb3RpZnMpW1syXV0KCnJvd25hbWVzKGRlZXBfbW90aWZfbXR4KSA8LSB0b2xvd2VyKHJvd25hbWVzKGRlZXBfbW90aWZfbXR4KSkKc3Vic3RyKHJvd25hbWVzKGRlZXBfbW90aWZfbXR4KSwgMSwgMSkgPC0gdG91cHBlcihzdWJzdHIocm93bmFtZXMoZGVlcF9tb3RpZl9tdHgpLCAxLCAxKSkKCgoKYGBgCgoKQ29sb3IgUGFsZXR0ZToKCmBgYHtyfQpjb2xQYWxldHRlX2NlbGx0eXBlcyA9IGMoJyM1MzJDOEEnLAogJyNjMTlmNzAnLAogJyNmOWRlY2YnLAogJyNjOWE5OTcnLAogJyNCNTFEOEQnLAogJyMzRjg0QUEnLAogJyM5ZTY3NjInLAogJyMzNTRFMjMnLAogJyNGMzk3QzAnLAogJyNmZjg5MWMnLAogJyM2MzU1NDcnLAogJyNDNzIyMjgnLAogJyNmNzkwODMnLAogJyNFRjRFMjInLAogJyM5ODk4OTgnLAogJyM3RjY4NzQnLAogJyM4ODcwYWQnLAogJyM2NDdhNGYnLAogJyNFRjVBOUQnLAogJyNGQkJFOTInLAogJyMxMzk5OTInLAogJyNjYzc4MTgnLAogJyNERkNERTQnLAogJyM4RUM3OTInLAogJyNDNTk0QkYnLAogJyNDM0MzODgnLAogJyMwRjRBOUMnLAogJyNGQUNCMTInLAogJyM4REI1Q0UnLAogJyMxQTFBMUEnLAogJyNDOUVCRkInLAogJyNEQUJFOTknLAogJyM2NUE4M0UnLAogJyMwMDU1NzknLAogJyNDREUwODgnLAogJyNmN2Y3OWUnLAogJyNGNkJGQ0InKQoKY2VsbHR5cGVzIDwtIChhcy5kYXRhLmZyYW1lKGdldENlbGxDb2xEYXRhKHByb2opKSAlPiUgZ3JvdXBfYnkoY2VsbHR5cGVzKSAlPiUgCiBzdW1tYXJpc2UobiA9IG4oKSkpJGNlbGx0eXBlcwoKY29sIDwtIHNldE5hbWVzKGNvbFBhbGV0dGVfY2VsbHR5cGVzLCBjZWxsdHlwZXMpCmBgYAoKClNlbGVjdCBhIGZldyBtYXJrZXIgZ2VuZXMgYW5kIHRyYW5zY3JpcHRpb24gZmFjdG9yczoKYGBge3J9Cm1hcmtlcl9nZW5lcyA8LSBjKCJMYW1iMSIsICAiU3BhcmMiLCAiRWxmNSIsICJBc2NsMiIsICJUZmFwMmMiLCAiVHRyIiwKICAgICAgICAgICAgICAgICAgIkFwb2EyIiwgIkFwb2UiLCAiQ3lzdG0xIiwgIkVtYiIsICJTcGluazEiLCAgIktydDE5IiwKICAgICAgICAgICAgICAgICAgIkRrazEiLCAiR3JobDMiLCAiVHJwNjMiLCAiR3JobDIiLCAgIlBheDYiLCAiUGF4MiIsCiAgICAgICAgICAgICAgICAgICJFbjEiLCAiRm94ZDMiLCAiVGZhcDJhIiwgIlBheDMiLCAiU294OSIsCiAgICAgICAgICAgICAgICAgICJTaXgzIiwgIkhlc3gxIiwgIklyeDMiLCAiSG94YjkiLCAiQ2R4NCIsCiAgICAgICAgICAgICAgICAgICJIZXMzIiwgIkhiYS1hMiIsICJIYmEtYTEiLCAgIkhiYi1iaDEiLCAiR2F0YTEiLAogICAgICAgICAgICAgICAgICAgIkdhdGE2IiwKICAgICAgICAgICAgICAgICAgIyJHYXRhNiIsICJHYXRhNSIsCiAgICAgICAgICAgICAgICAgICJDaXRlZDQiLAogICAgICAgICAgICAgICAgICAgIkNkaDUiLCAiUGVjYW0xIiwgIkFueGE1IiwgIkV0djIiLCAiSWdmMiIsCiAgICAgICAgICAgICAgICAgICJLcnQ4IiwgIktydDE4IiwgIlBtcDIyIiwgIkFobmFrIiwgIkJtcDQiLCAiVGJ4NCIsICJIb3hhMTEiLAogICAgICAgICAgICAgICAgICAiSG94YTEwIiwgIlRubnQyIiwgIk15bDQiLCAgIk15bDciLCAiQWN0YTIiLAogICAgICAgICAgICAgICAgICAiU21hcmNkMyIsICJUY2YyMSIsICJUYng2IiwgIkRsbDEiLCAiQWxkaDFhMiIsICJUY2YxNSIsCiAgICAgICAgICAgICAgICAgICJNZW94MSIsICJUYngxIiwgIkdieDIiLCAiQ2R4MSIsICJIb3hiMSIsICJIZXM3IiwgIk9zcjEiLAogICAgICAgICAgICAgICAgICAiTWVzcDIiLCAiTGVmdHkyIiwgIk1lc3AxIiwgIkNlcjEiLCAgIkNocmQiLCAKICAgICAgICAgICAgICAgICAgIkZveGEyIiwgIlBheDciLCAiRmdmOCIsICJMaHgxIiwgIk1peGwxIiwgIk90eDIiLCAiSGhleCIsCiAgICAgICAgICAgICAgICAgICAiSWZpdG0zIiwgIk5reDEtMiIsICJFb21lcyIsICJOYW5vZyIsICJVdGYxIiwKICAgICAgICAgICAgICAgICAgIkVwY2FtIiwgIlBvdTVmMSIgKQojIlNveDIiCiMiR2F0YTIiCiMiR2F0YTQiCiMgICJHYXRhNSIsCiMiR2F0YTYiLAojICJHc2MiLApgYGAKCmBgYCN7cn0KICBwIDwtIG1ldGFkYXRhICU+JQogICAgbXV0YXRlKCEhbiA6PSBtb3RpZl9uKSAlPiUKICAgIGdyb3VwX2J5KGNlbGx0eXBlcykgJT4lCiAgICAjc3VtbWFyaXNlX2F0KHZhcnMobiksIGZ1bnMobWVhbiguLCBuYS5ybT1UUlVFKSkpCiAgICBzdW1tYXJpc2UobWVhbiA9IG1lYW4oISEoc3ltKG4pKSkpICU+JQogICAgZ2dwbG90KCkgKwogICAgZ2VvbV9iYXIoYWVzKHggPSBjZWxsdHlwZXMsIHkgPSBtZWFuLCBmaWxsID0gY2VsbHR5cGVzKSwgc3RhdCA9ICJpZGVudGl0eSIpICsKICAgIHNjYWxlX2ZpbGxfbWFudWFsKHZhbHVlcyA9IGNvbCkgKwogICAgdGhlbWUoYXhpcy50ZXh0LnggPSBlbGVtZW50X3RleHQoYW5nbGUgPSA5MCwgdmp1c3QgPSAwLjUsIGhqdXN0PTEpLCBsZWdlbmQucG9zaXRpb24gPSAiTm9uZSIpICsjJT4lIHByaW50KCkKICAgIGxhYnModGl0bGUgPSBwYXN0ZTAobiksIHkgPSAiU0VBQ2VsbCBkZXZpYXRpb24gc2NvcmUiKQogIHByaW50KHApCmBgYAoKCiMgR0FUQSBGYWN0b3JzCgpgYGB7ciwgZmlnLndpZHRoPTEwLCBmaWcuaGVpZ2h0PTEwfQpmb3IgKG4gaW4gYygiR2F0YTEiLCAiR2F0YTIiLCAiR2F0YTMiLCAiR2F0YTQiLCAiR2F0YTYiKSkgeyAjICJHYXRhNSIsCiAgcHJpbnQobikKICAjIHNlbGVjdCBvbmUgcm93IGZvciBhIHBhcnRpY3VsYXIgZ2VuZQogIHNjb3JlX24gPC0gZ2VuZV9zY29yZXNfbWF0W3Jvd25hbWVzKGdlbmVfc2NvcmVzX21hdCkgJWluJSAgbiwgXQogICMgYWRkIHNjb3JlIGZvciB0aGlzIGdlbmUgdG8gbWV0YWRhdGEKICBtZXRhZGF0YVtwYXN0ZTAoInNjb3JlXyIsbildIDwtIHNjb3JlX24KICAjIHNlbGVjdCBnZW5lIGV4cHJlc3Npb24gZm9yIGEgcGFydGljdWxhciBnZW5lCiAgI2V4cHJfbiA8LSBsb2dub3JtW3Jvd25hbWVzKGxvZ25vcm0pICVpbiUgYyhuKSwgXQogICNtZXRhZGF0YVtwYXN0ZTAoImV4cHJfIiwgbildIDwtIGV4cHJfbgogIAogIAogIHNlYWNlbGxzX24gPC0gZGVlcF9tb3RpZl9tdHhbbiwgXQogIG1ldGFkYXRhW3Bhc3RlMCgic2VhY2VsbF8iLCBuKV0gPC0gc2VhY2VsbHNfbgoKICAjIHNlbGVjdCBtb3RpZiB6IHNjb3JlCiAgbW90aWZfbiA8LSBtb3RpZl9tdHhbcm93bmFtZXMobW90aWZfbXR4KSA9PSAgdGZzW2dyZXBsKHBhc3RlMCgiXiIsIG4pLCB0ZnMpXSwgXQogIG1ldGFkYXRhW3Bhc3RlMCgibW90aWZfIiwgbildIDwtIG1vdGlmX24KICAKICAjIGNyZWF0ZSBiYXJwbG90cyBmb3IgZ2VuZSBzY29yZXMsIGdlbmUgZXhwcmVzc2lvbiBhbmQgbW90aWYgei1zY29yZQogIHBsb3RzIDwtIG1hcChjKCJzY29yZV8iLCAibW90aWZfIiksIGZ1bmN0aW9uKHApewogICAgZGYgPC0gbWV0YWRhdGEgJT4lCiAgICBncm91cF9ieShjZWxsdHlwZXMpICU+JQogICAgc3VtbWFyaXNlX2F0KHZhcnMocGFzdGUwKHAsIG4pKSwgZnVucyhtZWFuKC4sIG5hLnJtPVRSVUUpKSkKICAgIGRmICU+JSBnZ3Bsb3QoKSArCiAgICBnZW9tX2JhcihhZXMoeCA9IGNlbGx0eXBlcywgeSA9IGRmICU+JSBwdWxsKHBhc3RlMChwLCBuKSksIAogICAgICAgICAgICAgICAgIGZpbGwgPSBjZWxsdHlwZXMpLCBzdGF0ID0gImlkZW50aXR5IikgKwogICAgc2NhbGVfZmlsbF9tYW51YWwodmFsdWVzID0gY29sKSArCiAgICB0aGVtZShheGlzLnRleHQueCA9IGVsZW1lbnRfdGV4dChhbmdsZSA9IDkwLCB2anVzdCA9IDAuNSwgaGp1c3Q9MSksCiAgICAgICAgICBsZWdlbmQucG9zaXRpb24gPSAibm9uZSIpICsKICAgIGxhYnModGl0bGUgPSBwYXN0ZTAobiksIHggPSAiY2VsbHR5cGUiLCB5ID0gcGFzdGUwKHApKQogIH0pCiAgc2VhX3Bsb3QgPC0gbWFwKHNlcS5pbnQoMSksIGZ1bmN0aW9uKHNlYSl7CiAgICBtZXRhZGF0YSAlPiUgCiAgICBncm91cF9ieShjZWxsdHlwZXMpICU+JQogICAgc3VtbWFyaXNlKG1lYW4gPSBtZWFuKCEhc3ltKHBhc3RlMCgic2VhY2VsbF8iLG4pKSkpICU+JQogICAgZ2dwbG90KCkgKwogICAgZ2VvbV9iYXIoYWVzKHggPSBjZWxsdHlwZXMsIHkgPSBtZWFuLCBmaWxsID0gY2VsbHR5cGVzKSwgc3RhdCA9ICJpZGVudGl0eSIpICsKICAgIHNjYWxlX2ZpbGxfbWFudWFsKHZhbHVlcyA9IGNvbCkgKwogICAgdGhlbWUoYXhpcy50ZXh0LnggPSBlbGVtZW50X3RleHQoYW5nbGUgPSA5MCwgdmp1c3QgPSAwLjUsIGhqdXN0PTEpLCBsZWdlbmQucG9zaXRpb24gPSAiTm9uZSIpICsjJT4lIHByaW50KCkKICAgIGxhYnModGl0bGUgPSBwYXN0ZTAobiksIHkgPSAiU0VBQ2VsbCBkZXZpYXRpb24gc2NvcmUiKQogIH0pCiAgIyBjcmVhdGUgc2NhdHRlciBwbG90cyBmb3IgZ2VuZSBleHByZXNzaW9uIGFuZCBtb3RpZiB6LXNjb3JlCiAgc2NhdHRlcl9wbG90cyA8LSBtYXAoc2VxLmludCgxKSwgZnVuY3Rpb24ocyl7CiAgICBkZiA8LSBtZXRhZGF0YSAlPiUgZ3JvdXBfYnkoY2VsbHR5cGVzKSAlPiUgCiAgICBzdW1tYXJpc2VfYXQodmFycyhwYXN0ZTAoIm1vdGlmXyIsIG4pLCAKICAgICAgICAgICAgICAgICAgICAgIHBhc3RlMCgic2NvcmVfIiwgbikpLCAKICAgICAgICAgICAgICAgICBmdW5zKG1lYW4oLiwgbmEucm0gPSBUUlVFKSkpCiAgICBkZiAlPiUKICAgIGdncGxvdCgpICsKICAgIGdlb21fc21vb3RoKGFlcyh4ID0gZGYgJT4lIHB1bGwocGFzdGUwKCJzY29yZV8iLCBuKSksIAogICAgICAgICAgICAgICAgICAgIHkgPSBkZiAlPiUgcHVsbChwYXN0ZTAoIm1vdGlmXyIsIG4pKSksCiAgICAgICAgICAgICAgICBmb3JtdWxhID0geSB+IHgsIG1ldGhvZCA9ICJsbSIsIHNpemUgPSAuMSkgKwogICAgZ2VvbV9wb2ludChhZXMoeCA9IGRmICU+JSBwdWxsKHBhc3RlMCgic2NvcmVfIiwgbikpLCAKICAgICAgICAgICAgICAgICAgIHkgPSBkZiAlPiUgcHVsbChwYXN0ZTAoIm1vdGlmXyIsIG4pKSwgCiAgICAgICAgICAgICAgICAgICBjb2wgPSBjZWxsdHlwZXMsIHNpemUgPSAxKSkgKwogICAgI2xhYnMoeCA9ICJHYXRhMSBnZW5lIGV4cHJlc3Npb24iLCB5ID0gIkdhdGExIG1vdGlmIGFjY2Vzc2liaWxpdHkgKHotc2NvcmUpIikgKwogICAgdGhlbWUobGVnZW5kLnBvc2l0aW9uID0gIm5vbmUiKSArCiAgICBzY2FsZV9jb2xvcl9tYW51YWwodmFsdWVzID0gY29sKSArCiAgICBsYWJzKHRpdGxlID0gcGFzdGUwKG4pLCB4ID0gImdlbmUgYWN0aXZpdHkgc2NvcmUiLCB5ID0gIlRGLW1vdGlmIHotc2NvcmUiKQogICAgfSkKICBkby5jYWxsKGdyaWRFeHRyYTo6Z3JpZC5hcnJhbmdlLCBjKHBsb3RzLCBzZWFfcGxvdCwgc2NhdHRlcl9wbG90cywgbmNvbCA9IDIpKQoKfQpgYGAKCiMgTWFya2VyIEdlbmVzCgojIyBTRUFDZWxscwoKCmBgYHtyLCBmaWcud2lkdGg9MTAsIGZpZy5oZWlnaHQ9MTB9Cgpmb3IgKG4gaW4gbWFya2VyX2dlbmVzKSB7CiAgcHJpbnQobikKICAjIHNlbGVjdCBvbmUgcm93IGZvciBhIHBhcnRpY3VsYXIgZ2VuZQogIHNjb3JlX24gPC0gZ2VuZV9zY29yZXNfbWF0W3Jvd25hbWVzKGdlbmVfc2NvcmVzX21hdCkgJWluJSBjKG4pLCBdCiAgIyBhZGQgc2NvcmUgZm9yIHRoaXMgZ2VuZSB0byBtZXRhZGF0YQogIG1ldGFkYXRhW3Bhc3RlMCgic2NvcmVfIixuKV0gPC0gc2NvcmVfbgogICMgc2VsZWN0IGdlbmUgZXhwcmVzc2lvbiBmb3IgYSBwYXJ0aWN1bGFyIGdlbmUKICBleHByX24gPC0gbG9nbm9ybVtyb3duYW1lcyhsb2dub3JtKSAlaW4lIGMobiksIF0KICBtZXRhZGF0YVtwYXN0ZTAoImV4cHJfIiwgbildIDwtIGV4cHJfbgogIAogIAogIHNlYWNlbGxzX24gPC0gc2VhX210eFtyb3duYW1lcyhzZWFfbXR4KSA9PSB0ZnNbZ3JlcGwocGFzdGUwKCJeIiwgbiksIHRmcyldLCBdCiAgc2VhX21ldGFbcGFzdGUwKCJzZWFjZWxsXyIsIG4pXSA8LSBzZWFjZWxsc19uCgogICMgaWYgdGhlIG1hcmtlciBnZW5lIGlzIGEgVEYKICBpZiAobGVuZ3RoKHRmc1tncmVwbChwYXN0ZTAoIl4iLCBuKSwgdGZzKV0pID4gMCkgewogICAgCiAgICAgICMgc2VsZWN0IG1vdGlmIHogc2NvcmUKICAgICAgbW90aWZfbiA8LSBtb3RpZl9tdHhbcm93bmFtZXMobW90aWZfbXR4KSA9PSAgdGZzW2dyZXBsKHBhc3RlMCgiXiIsIG4pLCB0ZnMpXSwgXQogICAgICBtZXRhZGF0YVtwYXN0ZTAoIm1vdGlmXyIsIG4pXSA8LSBtb3RpZl9uCiAgICAgIAogICAgICAjIGNyZWF0ZSBiYXJwbG90cyBmb3IgZ2VuZSBzY29yZXMsIGdlbmUgZXhwcmVzc2lvbiBhbmQgbW90aWYgei1zY29yZQogICAgICBwbG90cyA8LSBtYXAoYygic2NvcmVfIiwgIm1vdGlmXyIpLCBmdW5jdGlvbihwKXsKICAgICAgICBkZiA8LSBtZXRhZGF0YSAlPiUKICAgICAgICBncm91cF9ieShjZWxsdHlwZXMpICU+JQogICAgICAgIHN1bW1hcmlzZV9hdCh2YXJzKHBhc3RlMChwLCBuKSksIGZ1bnMobWVhbiguLCBuYS5ybT1UUlVFKSkpCiAgICAgICAgZGYgJT4lIGdncGxvdCgpICsKICAgICAgICBnZW9tX2JhcihhZXMoeCA9IGNlbGx0eXBlcywgeSA9IGRmICU+JSBwdWxsKHBhc3RlMChwLCBuKSksIAogICAgICAgICAgICAgICAgICAgICBmaWxsID0gY2VsbHR5cGVzKSwgc3RhdCA9ICJpZGVudGl0eSIpICsKICAgICAgICBzY2FsZV9maWxsX21hbnVhbCh2YWx1ZXMgPSBjb2wpICsKICAgICAgICB0aGVtZShheGlzLnRleHQueCA9IGVsZW1lbnRfdGV4dChhbmdsZSA9IDkwLCB2anVzdCA9IDAuNSwgaGp1c3Q9MSksCiAgICAgICAgICAgICAgbGVnZW5kLnBvc2l0aW9uID0gIm5vbmUiKSArCiAgICAgICAgbGFicyh0aXRsZSA9IHBhc3RlMChuKSwgeCA9ICJjZWxsdHlwZSIsIHkgPSBwYXN0ZTAocCkpCiAgICAgIH0pCiAgICAgIAogICAgICAjIGNyZWF0ZSBzY2F0dGVyIHBsb3RzIGZvciBnZW5lIGV4cHJlc3Npb24gYW5kIG1vdGlmIHotc2NvcmUKICAgICAgc2NhdHRlcl9wbG90cyA8LSBtYXAoc2VxLmludCgxKSwgZnVuY3Rpb24ocyl7CiAgICAgICAgZGYgPC0gbWV0YWRhdGEgJT4lIGdyb3VwX2J5KGNlbGx0eXBlcykgJT4lIAogICAgICAgIHN1bW1hcmlzZV9hdCh2YXJzKHBhc3RlMCgibW90aWZfIiwgbiksIAogICAgICAgICAgICAgICAgICAgICAgICAgIHBhc3RlMCgic2NvcmVfIiwgbikpLCAKICAgICAgICAgICAgICAgICAgICAgZnVucyhtZWFuKC4sIG5hLnJtID0gVFJVRSkpKQogICAgICAgIGRmICU+JQogICAgICAgIGdncGxvdCgpICsKICAgICAgICBnZW9tX3Ntb290aChhZXMoeCA9IGRmICU+JSBwdWxsKHBhc3RlMCgic2NvcmVfIiwgbikpLCAKICAgICAgICAgICAgICAgICAgICAgICAgeSA9IGRmICU+JSBwdWxsKHBhc3RlMCgibW90aWZfIiwgbikpKSwKICAgICAgICAgICAgICAgICAgICBmb3JtdWxhID0geSB+IHgsIG1ldGhvZCA9ICJsbSIsIHNpemUgPSAuMSkgKwogICAgICAgIGdlb21fcG9pbnQoYWVzKHggPSBkZiAlPiUgcHVsbChwYXN0ZTAoInNjb3JlXyIsIG4pKSwgCiAgICAgICAgICAgICAgICAgICAgICAgeSA9IGRmICU+JSBwdWxsKHBhc3RlMCgibW90aWZfIiwgbikpLCAKICAgICAgICAgICAgICAgICAgICAgICBjb2wgPSBjZWxsdHlwZXMsIHNpemUgPSAxKSkgKwogICAgICAgICNsYWJzKHggPSAiR2F0YTEgZ2VuZSBleHByZXNzaW9uIiwgeSA9ICJHYXRhMSBtb3RpZiBhY2Nlc3NpYmlsaXR5ICh6LXNjb3JlKSIpICsKICAgICAgICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAibm9uZSIpICsKICAgICAgICBzY2FsZV9jb2xvcl9tYW51YWwodmFsdWVzID0gY29sKSArCiAgICAgICAgbGFicyh0aXRsZSA9IHBhc3RlMChuKSwgeCA9ICJnZW5lIGFjdGl2aXR5IHNjb3JlIiwgeSA9ICJURi1tb3RpZiB6LXNjb3JlIikKICAgICAgICB9KQogICAgICAgIHNlYV9wbG90IDwtIG1hcChzZXEuaW50KDEpLCBmdW5jdGlvbihzZWEpewogICAgICBzZWFfbWV0YSAlPiUgCiAgICAgIGdyb3VwX2J5KGNlbGx0eXBlcykgJT4lCiAgICAgIHN1bW1hcmlzZShtZWFuID0gbWVhbighIXN5bShwYXN0ZTAoInNlYWNlbGxfIixuKSkpKSAlPiUKICAgICAgZ2dwbG90KCkgKwogICAgICBnZW9tX2JhcihhZXMoeCA9IGNlbGx0eXBlcywgeSA9IG1lYW4sIGZpbGwgPSBjZWxsdHlwZXMpLCBzdGF0ID0gImlkZW50aXR5IikgKwogICAgICBzY2FsZV9maWxsX21hbnVhbCh2YWx1ZXMgPSBjb2wpICsKICAgICAgdGhlbWUoYXhpcy50ZXh0LnggPSBlbGVtZW50X3RleHQoYW5nbGUgPSA5MCwgdmp1c3QgPSAwLjUsIGhqdXN0PTEpLCBsZWdlbmQucG9zaXRpb24gPSAiTm9uZSIpICsjJT4lIHByaW50KCkKICAgICAgbGFicyh0aXRsZSA9IHBhc3RlMChuKSwgeSA9ICJTRUFDZWxsIGRldmlhdGlvbiBzY29yZSIpCiAgICB9KQogICAgICAKICAgICMgaWYgdGhlIG1hcmtlciBnZW5lIGlzIG5vIFRGCiAgICB9IGVsc2UgewogICAgICAKICAgICAgIyBjcmVhdGUgYmFycGxvdHMgZm9yIGdlbmUgc2NvcmVzLCBnZW5lIGV4cHJlc3Npb24gYW5kIG1vdGlmIHotc2NvcmUKICAgICAgcGxvdHMgPC0gbWFwKGMoInNjb3JlXyIpLCBmdW5jdGlvbihwKXsKICAgICAgICBkZiA8LSBtZXRhZGF0YSAlPiUKICAgICAgICBncm91cF9ieShjZWxsdHlwZXMpICU+JQogICAgICAgIHN1bW1hcmlzZV9hdCh2YXJzKHBhc3RlMChwLCBuKSksIGZ1bnMobWVhbiguLCBuYS5ybT1UUlVFKSkpCiAgICAgICAgZGYgJT4lIGdncGxvdCgpICsKICAgICAgICBnZW9tX2JhcihhZXMoeCA9IGNlbGx0eXBlcywgeSA9IGRmICU+JSBwdWxsKHBhc3RlMChwLCBuKSksIAogICAgICAgICAgICAgICAgICAgICBmaWxsID0gY2VsbHR5cGVzKSwgc3RhdCA9ICJpZGVudGl0eSIpICsKICAgICAgICBzY2FsZV9maWxsX21hbnVhbCh2YWx1ZXMgPSBjb2wpICsKICAgICAgICB0aGVtZShheGlzLnRleHQueCA9IGVsZW1lbnRfdGV4dChhbmdsZSA9IDkwLCB2anVzdCA9IDAuNSwgaGp1c3Q9MSksCiAgICAgICAgICAgICAgbGVnZW5kLnBvc2l0aW9uID0gIm5vbmUiKSArCiAgICAgICAgbGFicyh0aXRsZSA9IHBhc3RlMChuKSwgeCA9ICJjZWxsdHlwZSIsIHkgPSBwYXN0ZTAocCkpCgogICAgICB9KQogICAgICAKICAgICAgIyBjcmVhdGUgc2NhdHRlciBwbG90cyBmb3IgZ2VuZSBleHByZXNzaW9uIGFuZCBnZW5lX3Njb3JlCiAgICAgIHNjYXR0ZXJfcGxvdHMgPC0gbWFwKHNlcS5pbnQoMSksIGZ1bmN0aW9uKHMpewogICAgICAgIGRmIDwtIG1ldGFkYXRhICU+JSBncm91cF9ieShjZWxsdHlwZXMpICU+JQogICAgICAgIHN1bW1hcmlzZV9hdCh2YXJzKHBhc3RlMCgic2NvcmVfIiwgbiksCiAgICAgICAgICAgICAgICAgICAgICAgICAgcGFzdGUwKCJleHByXyIsIG4pKSwKICAgICAgICAgICAgICAgICAgICAgZnVucyhtZWFuKC4sIG5hLnJtID0gVFJVRSkpKSAKICAgICAgICBkZiAlPiUKICAgICAgICBnZ3Bsb3QoKSArCiAgICAgICAgZ2VvbV9zbW9vdGgoYWVzKHggPSBkZiAlPiUgcHVsbChwYXN0ZTAoImV4cHJfIiwgbikpLCAKICAgICAgICAgICAgICAgICAgICAgICAgeSA9IGRmICU+JSBwdWxsKHBhc3RlMCgic2NvcmVfIiwgbikpKSwKICAgICAgICAgICAgICAgICAgICBmb3JtdWxhID0geSB+IHgsIG1ldGhvZCA9ICJsbSIsIHNpemUgPSAuMSkgKwogICAgICAgIGdlb21fcG9pbnQoYWVzKHggPSBkZiAlPiUgcHVsbChwYXN0ZTAoImV4cHJfIiwgbikpLAogICAgICAgICAgICAgICAgICAgICAgIHkgPSBkZiAlPiUgcHVsbChwYXN0ZTAoInNjb3JlXyIsIG4pKSwKICAgICAgICAgICAgICAgICAgICAgICBjb2wgPSBjZWxsdHlwZXMsIHNpemUgPSAxKSkgKwogICAgICAgICNsYWJzKHggPSAiR2F0YTEgZ2VuZSBleHByZXNzaW9uIiwgeSA9ICJHYXRhMSBtb3RpZiBhY2Nlc3NpYmlsaXR5ICh6LXNjb3JlKSIpICsKICAgICAgICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAibm9uZSIpICsKICAgICAgICBzY2FsZV9jb2xvcl9tYW51YWwodmFsdWVzID0gY29sKSsKICAgICAgICBsYWJzKHRpdGxlID0gcGFzdGUwKG4pLCB4ID0gImdlbmUgZXhwcmVzc2lvbiIsIHkgPSAiZ2VuZSBhY3Rpdml0eSBzY29yZSIpCgogICAgICB9KQogICAgfQoKICBkby5jYWxsKGdyaWRFeHRyYTo6Z3JpZC5hcnJhbmdlLCBjKHBsb3RzLCBzY2F0dGVyX3Bsb3RzLCBuY29sID0gMiwgbnJvdyA9IDIpKQogICMgCgp9CmBgYAoKCgoKUGxvdCBnZW5lIGV4cHJlc3Npb24sIGdlbmUgc2NvcmVzICYgbW90aWYgei1zY29yZXM6CgpgYGAje3IsIGZpZy53aWR0aD0xMCwgZmlnLmhlaWdodD0xMH0KCmZvciAobiBpbiBtYXJrZXJfZ2VuZXMpIHsKICBwcmludChuKQogICMgc2VsZWN0IG9uZSByb3cgZm9yIGEgcGFydGljdWxhciBnZW5lCiAgc2NvcmVfbiA8LSBnZW5lX3Njb3Jlc19tYXRbcm93bmFtZXMoZ2VuZV9zY29yZXNfbWF0KSAlaW4lIGMobiksIF0KICAjIGFkZCBzY29yZSBmb3IgdGhpcyBnZW5lIHRvIG1ldGFkYXRhCiAgbWV0YWRhdGFbcGFzdGUwKCJzY29yZV8iLG4pXSA8LSBzY29yZV9uCiAgIyBzZWxlY3QgZ2VuZSBleHByZXNzaW9uIGZvciBhIHBhcnRpY3VsYXIgZ2VuZQogIGV4cHJfbiA8LSBsb2dub3JtW3Jvd25hbWVzKGxvZ25vcm0pICVpbiUgYyhuKSwgXQogIG1ldGFkYXRhW3Bhc3RlMCgiZXhwcl8iLCBuKV0gPC0gZXhwcl9uCiAgCiAgCiAgc2VhY2VsbHNfbiA8LSBzZWFfbXR4W3Jvd25hbWVzKHNlYV9tdHgpID09IHRmc1tncmVwbChwYXN0ZTAoIl4iLCBuKSwgdGZzKV0sIF0KICBzZWFfbWV0YVtwYXN0ZTAoInNlYWNlbGxfIiwgbildIDwtIHNlYWNlbGxzX24KCiAgIyBpZiB0aGUgbWFya2VyIGdlbmUgaXMgYSBURgogIGlmIChsZW5ndGgodGZzW2dyZXBsKHBhc3RlMCgiXiIsIG4pLCB0ZnMpXSkgPiAwKSB7CiAgICAKICAgICAgIyBzZWxlY3QgbW90aWYgeiBzY29yZQogICAgICBtb3RpZl9uIDwtIG1vdGlmX210eFtyb3duYW1lcyhtb3RpZl9tdHgpID09ICB0ZnNbZ3JlcGwocGFzdGUwKCJeIiwgbiksIHRmcyldLCBdCiAgICAgIG1ldGFkYXRhW3Bhc3RlMCgibW90aWZfIiwgbildIDwtIG1vdGlmX24KICAgICAgCiAgICAgICMgY3JlYXRlIGJhcnBsb3RzIGZvciBnZW5lIHNjb3JlcywgZ2VuZSBleHByZXNzaW9uIGFuZCBtb3RpZiB6LXNjb3JlCiAgICAgIHBsb3RzIDwtIG1hcChjKCJzY29yZV8iLCAiZXhwcl8iLCAibW90aWZfIiksIGZ1bmN0aW9uKHApewogICAgICAgIGRmIDwtIG1ldGFkYXRhICU+JQogICAgICAgIGdyb3VwX2J5KGNlbGx0eXBlcykgJT4lCiAgICAgICAgc3VtbWFyaXNlX2F0KHZhcnMocGFzdGUwKHAsIG4pKSwgZnVucyhtZWFuKC4sIG5hLnJtPVRSVUUpKSkKICAgICAgICBkZiAlPiUgZ2dwbG90KCkgKwogICAgICAgIGdlb21fYmFyKGFlcyh4ID0gY2VsbHR5cGVzLCB5ID0gZGYgJT4lIHB1bGwocGFzdGUwKHAsIG4pKSwgCiAgICAgICAgICAgICAgICAgICAgIGZpbGwgPSBjZWxsdHlwZXMpLCBzdGF0ID0gImlkZW50aXR5IikgKwogICAgICAgIHNjYWxlX2ZpbGxfbWFudWFsKHZhbHVlcyA9IGNvbCkgKwogICAgICAgIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KGFuZ2xlID0gOTAsIHZqdXN0ID0gMC41LCBoanVzdD0xKSwKICAgICAgICAgICAgICBsZWdlbmQucG9zaXRpb24gPSAibm9uZSIpICsKICAgICAgICBsYWJzKHRpdGxlID0gcGFzdGUwKG4pLCB4ID0gImNlbGx0eXBlIiwgeSA9IHBhc3RlMChwKSkKICAgICAgfSkKICAgICAgCiAgICAgICMgY3JlYXRlIHNjYXR0ZXIgcGxvdHMgZm9yIGdlbmUgZXhwcmVzc2lvbiBhbmQgbW90aWYgei1zY29yZQogICAgICBzY2F0dGVyX3Bsb3RzIDwtIG1hcChzZXEuaW50KDEpLCBmdW5jdGlvbihzKXsKICAgICAgICBkZiA8LSBtZXRhZGF0YSAlPiUgZ3JvdXBfYnkoY2VsbHR5cGVzKSAlPiUgCiAgICAgICAgc3VtbWFyaXNlX2F0KHZhcnMocGFzdGUwKCJtb3RpZl8iLCBuKSwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgcGFzdGUwKCJzY29yZV8iLCBuKSksIAogICAgICAgICAgICAgICAgICAgICBmdW5zKG1lYW4oLiwgbmEucm0gPSBUUlVFKSkpCiAgICAgICAgZGYgJT4lCiAgICAgICAgZ2dwbG90KCkgKwogICAgICAgIGdlb21fc21vb3RoKGFlcyh4ID0gZGYgJT4lIHB1bGwocGFzdGUwKCJzY29yZV8iLCBuKSksIAogICAgICAgICAgICAgICAgICAgICAgICB5ID0gZGYgJT4lIHB1bGwocGFzdGUwKCJtb3RpZl8iLCBuKSkpLAogICAgICAgICAgICAgICAgICAgIGZvcm11bGEgPSB5IH4geCwgbWV0aG9kID0gImxtIiwgc2l6ZSA9IC4xKSArCiAgICAgICAgZ2VvbV9wb2ludChhZXMoeCA9IGRmICU+JSBwdWxsKHBhc3RlMCgic2NvcmVfIiwgbikpLCAKICAgICAgICAgICAgICAgICAgICAgICB5ID0gZGYgJT4lIHB1bGwocGFzdGUwKCJtb3RpZl8iLCBuKSksIAogICAgICAgICAgICAgICAgICAgICAgIGNvbCA9IGNlbGx0eXBlcywgc2l6ZSA9IDEpKSArCiAgICAgICAgI2xhYnMoeCA9ICJHYXRhMSBnZW5lIGV4cHJlc3Npb24iLCB5ID0gIkdhdGExIG1vdGlmIGFjY2Vzc2liaWxpdHkgKHotc2NvcmUpIikgKwogICAgICAgIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICJub25lIikgKwogICAgICAgIHNjYWxlX2NvbG9yX21hbnVhbCh2YWx1ZXMgPSBjb2wpICsKICAgICAgICBsYWJzKHRpdGxlID0gcGFzdGUwKG4pLCB4ID0gImdlbmUgYWN0aXZpdHkgc2NvcmUiLCB5ID0gIlRGLW1vdGlmIHotc2NvcmUiKQogICAgICAgIH0pCiAgICAgIAogICAgIyBpZiB0aGUgbWFya2VyIGdlbmUgaXMgbm8gVEYKICAgIH0gZWxzZSB7CiAgICAgIAogICAgICAjIGNyZWF0ZSBiYXJwbG90cyBmb3IgZ2VuZSBzY29yZXMsIGdlbmUgZXhwcmVzc2lvbiBhbmQgbW90aWYgei1zY29yZQogICAgICBwbG90cyA8LSBtYXAoYygic2NvcmVfIiwgImV4cHJfIiksIGZ1bmN0aW9uKHApewogICAgICAgIGRmIDwtIG1ldGFkYXRhICU+JQogICAgICAgIGdyb3VwX2J5KGNlbGx0eXBlcykgJT4lCiAgICAgICAgc3VtbWFyaXNlX2F0KHZhcnMocGFzdGUwKHAsIG4pKSwgZnVucyhtZWFuKC4sIG5hLnJtPVRSVUUpKSkKICAgICAgICBkZiAlPiUgZ2dwbG90KCkgKwogICAgICAgIGdlb21fYmFyKGFlcyh4ID0gY2VsbHR5cGVzLCB5ID0gZGYgJT4lIHB1bGwocGFzdGUwKHAsIG4pKSwgCiAgICAgICAgICAgICAgICAgICAgIGZpbGwgPSBjZWxsdHlwZXMpLCBzdGF0ID0gImlkZW50aXR5IikgKwogICAgICAgIHNjYWxlX2ZpbGxfbWFudWFsKHZhbHVlcyA9IGNvbCkgKwogICAgICAgIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KGFuZ2xlID0gOTAsIHZqdXN0ID0gMC41LCBoanVzdD0xKSwKICAgICAgICAgICAgICBsZWdlbmQucG9zaXRpb24gPSAibm9uZSIpICsKICAgICAgICBsYWJzKHRpdGxlID0gcGFzdGUwKG4pLCB4ID0gImNlbGx0eXBlIiwgeSA9IHBhc3RlMChwKSkKCiAgICAgIH0pCiAgICAgIAogICAgICAjIGNyZWF0ZSBzY2F0dGVyIHBsb3RzIGZvciBnZW5lIGV4cHJlc3Npb24gYW5kIGdlbmVfc2NvcmUKICAgICAgc2NhdHRlcl9wbG90cyA8LSBtYXAoc2VxLmludCgxKSwgZnVuY3Rpb24ocyl7CiAgICAgICAgZGYgPC0gbWV0YWRhdGEgJT4lIGdyb3VwX2J5KGNlbGx0eXBlcykgJT4lCiAgICAgICAgc3VtbWFyaXNlX2F0KHZhcnMocGFzdGUwKCJzY29yZV8iLCBuKSwKICAgICAgICAgICAgICAgICAgICAgICAgICBwYXN0ZTAoImV4cHJfIiwgbikpLAogICAgICAgICAgICAgICAgICAgICBmdW5zKG1lYW4oLiwgbmEucm0gPSBUUlVFKSkpIAogICAgICAgIGRmICU+JQogICAgICAgIGdncGxvdCgpICsKICAgICAgICBnZW9tX3Ntb290aChhZXMoeCA9IGRmICU+JSBwdWxsKHBhc3RlMCgiZXhwcl8iLCBuKSksIAogICAgICAgICAgICAgICAgICAgICAgICB5ID0gZGYgJT4lIHB1bGwocGFzdGUwKCJzY29yZV8iLCBuKSkpLAogICAgICAgICAgICAgICAgICAgIGZvcm11bGEgPSB5IH4geCwgbWV0aG9kID0gImxtIiwgc2l6ZSA9IC4xKSArCiAgICAgICAgZ2VvbV9wb2ludChhZXMoeCA9IGRmICU+JSBwdWxsKHBhc3RlMCgiZXhwcl8iLCBuKSksCiAgICAgICAgICAgICAgICAgICAgICAgeSA9IGRmICU+JSBwdWxsKHBhc3RlMCgic2NvcmVfIiwgbikpLAogICAgICAgICAgICAgICAgICAgICAgIGNvbCA9IGNlbGx0eXBlcywgc2l6ZSA9IDEpKSArCiAgICAgICAgI2xhYnMoeCA9ICJHYXRhMSBnZW5lIGV4cHJlc3Npb24iLCB5ID0gIkdhdGExIG1vdGlmIGFjY2Vzc2liaWxpdHkgKHotc2NvcmUpIikgKwogICAgICAgIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICJub25lIikgKwogICAgICAgIHNjYWxlX2NvbG9yX21hbnVhbCh2YWx1ZXMgPSBjb2wpKwogICAgICAgIGxhYnModGl0bGUgPSBwYXN0ZTAobiksIHggPSAiZ2VuZSBleHByZXNzaW9uIiwgeSA9ICJnZW5lIGFjdGl2aXR5IHNjb3JlIikKCiAgICAgIH0pCiAgICB9CgogIGRvLmNhbGwoZ3JpZEV4dHJhOjpncmlkLmFycmFuZ2UsIGMocGxvdHMsIHNjYXR0ZXJfcGxvdHMsIG5jb2wgPSAyLCBucm93ID0gMikpCiAgIyAKCn0KYGBgCgoKCmBgYCN7ciwgZmlnLndpZHRoPTEwLCBmaWcuaGVpZ2h0PTEwfQpwMSA8LSBtZXRhZGF0YSAlPiUKICBncm91cF9ieShjZWxsdHlwZXMpICU+JQogIHN1bW1hcmlzZV9hdCh2YXJzKEdhdGExKSwgZnVucyhtZWRpYW4oLiwgbmEucm09VFJVRSkpKSAlPiUgZ2dwbG90KCkgKwogIGdlb21fYmFyKGFlcyh4ID0gY2VsbHR5cGVzLCB5ID0gR2F0YTEsIGZpbGwgPSBjZWxsdHlwZXMpLCBzdGF0ID0gImlkZW50aXR5IikgKwogIHNjYWxlX2ZpbGxfbWFudWFsKHZhbHVlcyA9IGNvbCkgKwogIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KGFuZ2xlID0gOTAsIHZqdXN0ID0gMC41LCBoanVzdD0xKSwgCiAgICAgICAgbGVnZW5kLnBvc2l0aW9uID0gIm5vbmUiKQoKcDIgPC0gbWV0YWRhdGEgJT4lCiAgZ3JvdXBfYnkoY2VsbHR5cGVzKSAlPiUKICBzdW1tYXJpc2VfYXQodmFycyhnYXRhMV9zY29yZV9wMmcpLCBmdW5zKG1lYW4oLiwgbmEucm09VFJVRSkpKSAlPiUgZ2dwbG90KCkgKwogIGdlb21fYmFyKGFlcyh4ID0gY2VsbHR5cGVzLCB5ID0gZ2F0YTFfc2NvcmVfcDJnLCBmaWxsID0gY2VsbHR5cGVzKSwgc3RhdCA9ICJpZGVudGl0eSIpICsKICBzY2FsZV9maWxsX21hbnVhbCh2YWx1ZXMgPSBjb2wpICsKICB0aGVtZShheGlzLnRleHQueCA9IGVsZW1lbnRfdGV4dChhbmdsZSA9IDkwLCB2anVzdCA9IDAuNSwgaGp1c3Q9MSksIAogICAgICAgIGxlZ2VuZC5wb3NpdGlvbiA9ICJub25lIikKCgpwMyA8LSBtZXRhZGF0YSAlPiUgZ3JvdXBfYnkoY2VsbHR5cGVzKSAlPiUgCiAgc3VtbWFyaXNlX2F0KHZhcnMoZ2F0YTFfel9zY29yZXMsIEdhdGExKSwgZnVucyhtZWFuKC4sIG5hLnJtID0gVFJVRSkpKSAlPiUKICBnZ3Bsb3QoKSArCiAgZ2VvbV9zbW9vdGgoYWVzKHggPSBHYXRhMSwgeSA9IGdhdGExX3pfc2NvcmVzKSwgZm9ybXVsYSA9IHkgfiB4LCBtZXRob2QgPSAibG0iLCBzaXplID0gLjEpICsKICBnZW9tX3BvaW50KGFlcyh4ID0gR2F0YTEsIHkgPSBnYXRhMV96X3Njb3JlcywgY29sID0gY2VsbHR5cGVzLCBzaXplID0gMSkpICsKICBsYWJzKHggPSAiR2F0YTEgZ2VuZSBleHByZXNzaW9uIiwgeSA9ICJHYXRhMSBtb3RpZiBhY2Nlc3NpYmlsaXR5ICh6LXNjb3JlKSIpICsKICB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAibm9uZSIpICsKICBzY2FsZV9jb2xvcl9tYW51YWwodmFsdWVzID0gY29sKSAKCnA0IDwtIG1ldGFkYXRhICU+JQogIGdyb3VwX2J5KGNlbGx0eXBlcykgJT4lCiAgc3VtbWFyaXNlX2F0KHZhcnMoZ2F0YTFfel9zY29yZXMpLCBmdW5zKG1lYW4oLiwgbmEucm09VFJVRSkpKSAlPiUgZ2dwbG90KCkgKwogIGdlb21fYmFyKGFlcyh4ID0gY2VsbHR5cGVzLCB5ID0gZ2F0YTFfel9zY29yZXMsIGZpbGwgPSBjZWxsdHlwZXMpLCBzdGF0ID0gImlkZW50aXR5IikgKwogIHNjYWxlX2ZpbGxfbWFudWFsKHZhbHVlcyA9IGNvbCkgKwogIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KGFuZ2xlID0gOTAsIHZqdXN0ID0gMC41LCBoanVzdD0xKSwgCiAgICAgICAgbGVnZW5kLnBvc2l0aW9uID0gIm5vbmUiKQoKY293cGxvdDo6cGxvdF9ncmlkKHAxLCBwMiwgcDMsIHA0KQoKcDEgPC0gbWV0YWRhdGEgJT4lCiAgZ3JvdXBfYnkoY2VsbHR5cGVzKSAlPiUKICBzdW1tYXJpc2VfYXQodmFycyhTb3g5KSwgZnVucyhtZWFuKC4sIG5hLnJtPVRSVUUpKSkgJT4lIGdncGxvdCgpICsKICBnZW9tX2JhcihhZXMoeCA9IGNlbGx0eXBlcywgeSA9IFNveDksIGZpbGwgPSBjZWxsdHlwZXMpLCBzdGF0ID0gImlkZW50aXR5IikgKwogIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KGFuZ2xlID0gOTAsIHZqdXN0ID0gMC41LCBoanVzdD0xKSwKICAgICAgICBsZWdlbmQucG9zaXRpb24gPSAibm9uZSIpICsKICBzY2FsZV9maWxsX21hbnVhbCh2YWx1ZXMgPSBjb2wpCgoKCnAyIDwtIG1ldGFkYXRhICU+JQogIGdyb3VwX2J5KGNlbGx0eXBlcykgJT4lCiAgc3VtbWFyaXNlX2F0KHZhcnMoc294OV96X3Njb3JlcyksIGZ1bnMobWVhbiguLCBuYS5ybT1UUlVFKSkpICU+JSBnZ3Bsb3QoKSArCiAgZ2VvbV9iYXIoYWVzKHggPSBjZWxsdHlwZXMsIHkgPSBzb3g5X3pfc2NvcmVzLCBmaWxsID0gY2VsbHR5cGVzKSwgc3RhdCA9ICJpZGVudGl0eSIpICsKICBzY2FsZV9maWxsX21hbnVhbCh2YWx1ZXMgPSBjb2wpICsKICB0aGVtZShheGlzLnRleHQueCA9IGVsZW1lbnRfdGV4dChhbmdsZSA9IDkwLCB2anVzdCA9IDAuNSwgaGp1c3Q9MSksIAogICAgICAgIGxlZ2VuZC5wb3NpdGlvbiA9ICJub25lIikKCnAzIDwtIG1ldGFkYXRhICU+JQogIGdyb3VwX2J5KGNlbGx0eXBlcykgJT4lCiAgc3VtbWFyaXNlX2F0KHZhcnMoc294OV9zY29yZV9wMmcpLCBmdW5zKG1lYW4oLiwgbmEucm09VFJVRSkpKSAlPiUgZ2dwbG90KCkgKwogIGdlb21fYmFyKGFlcyh4ID0gY2VsbHR5cGVzLCB5ID0gc294OV9zY29yZV9wMmcsIGZpbGwgPSBjZWxsdHlwZXMpLCBzdGF0ID0gImlkZW50aXR5IikgKwogIHNjYWxlX2ZpbGxfbWFudWFsKHZhbHVlcyA9IGNvbCkgKwogIHRoZW1lKGF4aXMudGV4dC54ID0gZWxlbWVudF90ZXh0KGFuZ2xlID0gOTAsIHZqdXN0ID0gMC41LCBoanVzdD0xKSwgCiAgICAgICAgbGVnZW5kLnBvc2l0aW9uID0gIm5vbmUiKQoKCnA0IDwtIG1ldGFkYXRhICU+JSBncm91cF9ieShjZWxsdHlwZXMpICU+JSAKICBzdW1tYXJpc2VfYXQodmFycyhzb3g5X3pfc2NvcmVzLCBTb3g5KSwgZnVucyhtZWFuKC4sIG5hLnJtID0gVFJVRSkpKSAlPiUKICBnZ3Bsb3QoKSArCiAgZ2VvbV9zbW9vdGgoYWVzKHggPSBTb3g5LCB5ID0gc294OV96X3Njb3JlcyksIGZvcm11bGEgPSB5IH4geCwgbWV0aG9kID0gImxtIiwgc2l6ZSA9IC4xKSArCiAgZ2VvbV9wb2ludChhZXMoeCA9IFNveDksIHkgPSBzb3g5X3pfc2NvcmVzLCBjb2wgPSBjZWxsdHlwZXMsIHNpemUgPSAxKSkgKwogIGxhYnMoeCA9ICJHYXRhMSBnZW5lIGV4cHJlc3Npb24iLCB5ID0gIkdhdGExIG1vdGlmIGFjY2Vzc2liaWxpdHkgKHotc2NvcmUpIikgKwogIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICJub25lIikgKwogIHNjYWxlX2NvbG9yX21hbnVhbCh2YWx1ZXMgPSBjb2wpIAoKY293cGxvdDo6cGxvdF9ncmlkKHAxLCBwMiwgcDMsIHA0KQoKYGBgCgoKCgoKCg==